Introduction
First get the raw data into R and tidy it.
Reproducing figure 1b through 1g
Data transformation
Spectral analyses
Reproducing Table 1 using ELE supplement data.
Predictability (Figure 2)
Lyapunov exponents by direct method (Figure 3)
Lyapunov exponents by indirect method

Introduction

This is an attempt to reproduce the anaylses presented in the paper Chaos in a long-term experiment with a plankton community, by Elisa Benincà and others (the paper on the Nature website). Details of the methods are in the Supplement to the Nature paper.

Here include some text about what is presented below in the basic reproduction, and what additional was done and where that is described.

The data are available as an Excel file supplement to an Ecology Letters publication. The Excel file contains several datasheets. Two are particularly important, as they are the source of the raw data (one contains original species abundances, the one with the nutrient concentrations). Another datasheet in the ELE supplement contains transformed variables.

First get the raw data into R and tidy it.

rm(list=ls())
library(tidyr)
library(dplyr)

## 
## Attaching package: 'dplyr'
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(lubridate)
library(stringr)
library(ggplot2)
library(RCurl)

## Loading required package: bitops

library(pracma)
library(oce)

## Loading required package: mapproj
## Loading required package: maps
## 
## Attaching package: 'oce'
## 
## The following objects are masked from 'package:pracma':
## 
##     detrend, grad
## 
## The following object is masked from 'package:tidyr':
## 
##     extract

spp.abund <- read.csv(text=getURL("https://raw.githubusercontent.com/opetchey/RREEBES/master/Beninca_etal_2008_Nature/data/species_abundances_original.csv"), skip=7, header=T)

spp.abund <- select(spp.abund, -X, -X.1)
spp.abund <- spp.abund[-804:-920,]
str(spp.abund)

## 'data.frame':    803 obs. of  12 variables:
##  $ Date               : Factor w/ 803 levels "","01/02/91",..: 306 440 465 498 520 601 628 673 699 778 ...
##  $ Day.number         : int  1 6 7 8 9 12 13 15 16 19 ...
##  $ Cyclopoids         : num  0 0 0.0353 0 0.0353 ...
##  $ Calanoid.copepods  : num  1.04 2.03 1.72 2.41 1.71 ...
##  $ Rotifers           : num  7.7 10.19 8.08 6.06 5.94 ...
##  $ Protozoa           : Factor w/ 330 levels "","0","0,000001",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Nanophytoplankton  : num  0.106 0.212 0.212 0.212 0.212 ...
##  $ Picophytoplankton  : num  1 2 1.52 1.52 1.98 ...
##  $ Filamentous.diatoms: num  0 0 0 0 0 0 0 0 0 0 ...
##  $ Ostracods          : num  0 0 0 0.0187 0 ...
##  $ Harpacticoids      : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ Bacteria           : num  2.15 1.97 1.79 1.61 1.43 ...

The Protozoa variable contains some numbers with comman as the decimal separator. This creates a question about what dataset was used for the original analyses, as it could not have been this one.

spp.abund$Protozoa <- as.numeric(str_replace(spp.abund$Protozoa, ",", "."))

Format the dates as dates

spp.abund$Date <- dmy(spp.abund$Date)

Ooops… R assumes the experiment was done in the 21st century. Shouldn’t matter too much.

Check dates match the Day.number (should give true):

sum(spp.abund$Day.number == 1+as.numeric((spp.abund$Date - spp.abund$Date[1]) / 24 / 60 / 60)) == length(spp.abund$Date)

## [1] TRUE

Check for duplicate dates:

spp.abund$Date[duplicated(spp.abund$Date)]

## [1] "2096-10-28 UTC"

which(duplicated(spp.abund$Date))

## [1] 702

Original dataset contains a duplicated date: 28/10/1996 (row 709 and 710 in excel sheet). Lets change the date in row 709 to 26/10/1996, which will put it half way between the two surrounding dates:

which(spp.abund$Date==ymd("2096-10-28 UTC"))

## [1] 701 702

spp.abund$Date[701] <- ymd("2096-10-26 UTC")

Check dates match the Day.number (should give true):

sum(spp.abund$Day.number == 1+as.numeric((spp.abund$Date - spp.abund$Date[1]) / 24 / 60 / 60)) == length(spp.abund$Date)

## [1] FALSE

Fix the Day.number problem:

spp.abund$Day.number <- 1+as.numeric((spp.abund$Date - spp.abund$Date[1]) / 24 / 60 / 60)

Data is in wide format, so change it to long:

spp.abund <- gather(spp.abund, "variable", "value", 3:12)
str(spp.abund)

## 'data.frame':    8030 obs. of  4 variables:
##  $ Date      : POSIXct, format: "2090-07-12" "2090-07-17" ...
##  $ Day.number: num  1 6 7 8 9 12 13 15 16 19 ...
##  $ variable  : Factor w/ 10 levels "Cyclopoids","Calanoid.copepods",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ value     : num  0 0 0.0353 0 0.0353 ...

Bring in the nutrient data:

nuts <- read.csv(text=getURL("https://raw.githubusercontent.com/opetchey/RREEBES/master/Beninca_etal_2008_Nature/data/nutrients_original.csv"), skip=7, header=T)
nuts <- select(nuts, -X, -X.1)
nuts <- nuts[-349:-8163,]
nuts$Date <- dmy(nuts$Date)
nuts <- select(nuts, -NO2, -NO3, -NH4)
nuts$Date[duplicated(nuts$Date)]

## character(0)

which(duplicated(nuts$Date))

## integer(0)

nuts <- gather(nuts, "variable", "value", 3:4)
str(nuts)

## 'data.frame':    696 obs. of  4 variables:
##  $ Date      : POSIXct, format: "2090-09-18" "2090-09-24" ...
##  $ Day.number: int  69 75 82 90 96 103 110 117 124 131 ...
##  $ variable  : Factor w/ 2 levels "Total.dissolved.inorganic.nitrogen",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ value     : num  28.32 20.84 11.15 15.5 5.92 ...

Now put the two datasets together

all.data <- rbind(spp.abund, nuts)

Now select only the date range used in the Nature paper. From the supplment The analysis in Benincà et al. (Nature 2008) covered all data from 16/06/1991 until 20/10/1997. (Remembering dates in the R dataframes are 2090s.)

all.data <- filter(all.data, Date>dmy("16/06/2091") & Date<dmy("20/10/2097"))

Reproducing figure 1b through 1g

(No attempt to reproduce Figure 1a, as its a food web diagram.)

First quick go:

ggplot(all.data, aes(x=Day.number, y=value)) + geom_line() +
  facet_wrap(~variable, scales="free_y")

The code from here on needs modifying, as the object names are Owen’s old ones, and need changing to those used above.

Now we add a column that gives the variable types, same as in figure 1b through 1g. First make a lookup table giving species type:

tt <- data.frame(variable=unique(all.data$variable),
                 Type=c("Cyclopoids", "Herbivore", "Herbivore", "Herbivore",
                        "Phytoplankton",  "Phytoplankton", "Phytoplankton",
                        "Detritivore", "Detritivore", "Bacteria", "Nutrient", "Nutrient"))
tt

##                              variable          Type
## 1                          Cyclopoids    Cyclopoids
## 2                   Calanoid.copepods     Herbivore
## 3                            Rotifers     Herbivore
## 4                            Protozoa     Herbivore
## 5                   Nanophytoplankton Phytoplankton
## 6                   Picophytoplankton Phytoplankton
## 7                 Filamentous.diatoms Phytoplankton
## 8                           Ostracods   Detritivore
## 9                       Harpacticoids   Detritivore
## 10                           Bacteria      Bacteria
## 11 Total.dissolved.inorganic.nitrogen      Nutrient
## 12        Soluble.reactive.phosphorus      Nutrient

And add the Type variable to the new dataset:

all.data <- merge(all.data, tt)

First lets set the colours as in the original:

species.colour.mapping <- c("Cyclopoids"="pink",
                            "Calanoid.copepods"="red",
                            "Rotifers"="blue",
                            "Protozoa"="green",
                            "Nanophytoplankton"="red",
                            "Picophytoplankton"="black",
                            "Filamentous.diatoms"="green",
                            "Ostracods"="lightblue",
                            "Harpacticoids"="purple",
                            "Bacteria"="black",
                            "Total.dissolved.inorganic.nitrogen"="red",
                            "Soluble.reactive.phosphorus"="black")

Next change the order of the levels in the Type variable, so plots appear in the same order as in the original figure:

all.data$Type <- factor(all.data$Type, levels=c("Cyclopoids", "Herbivore", "Phytoplankton", "Nutrient",
                                    "Bacteria", "Detritivore"))

Now a version that doesn’t try to recreate the “gap” in the y axes of the original figures:

g1 <- qplot(as.numeric(Day.number), value, col=variable, data=all.data) +
  facet_wrap(~Type, ncol=2, scales="free_y") +
  geom_point() + geom_line() +
  scale_colour_manual(values = species.colour.mapping)
g1

Looks reasonably good.

Now a version that approximates the “gap”, by removing data above it:

an2 <- filter(all.data, Type=="Cyclopoids" & value<0.6 |
                Type=="Herbivore" & value<13 |
                Type=="Phytoplankton" & value<1400 |
                Type=="Nutrient" & value<50 |
                Type=="Bacteria" & value<10 |
                Type=="Detritivore" & value<0.7) 
g1 <- qplot(as.numeric(Day.number), value, col=variable, data=an2) +
  facet_wrap(~Type, ncol=2, scales="free_y") +
  geom_point() + geom_line() +
  scale_colour_manual(values = species.colour.mapping)
g1

Difficult it look like the data go off the top of the graph in ggplot.

Try logarithmic y-axes:

g1 <- qplot(as.numeric(Day.number), log10(value+0.00001), col=variable, data=all.data) +
  facet_wrap(~Type, ncol=2, scales="free_y") +
  geom_point() + geom_line() +
  scale_colour_manual(values = species.colour.mapping)
g1

Try fourth root, as this is used in the subsequent stats.

g1 <- qplot(as.numeric(Day.number), value^0.25, col=variable, data=all.data) +
  facet_wrap(~Type, ncol=2, scales="free_y") +
  geom_point() + geom_line() +
  scale_colour_manual(values = species.colour.mapping)
g1

Data transformation

Now we need to work with transformed data. Details of the transformation, copied from the Supplmentary information are in indented quote style in the following sections… looks like this:

Transformation of the time series. We transformed the original time series, shown in Fig. 1b-g of the main text, to obtain stationary time series with equidistant data and homogeneous units of measurement. The transformation steps are illustrated for the bacteria (Fig. S1).

Aside: The ELE supplement contains the raw data and the transformed data, in separate data sheets. I (Owen) also got the interpolated data from Stephen Ellner directly.

Interpolation

First, the time series were interpolated using cubic hermite interpolation, to obtain data with equidistant time intervals of 3.35 days (Fig. S1a).

Make a sequence of times at which to interpolate.

#aggregate(Day.number ~ variable, all.data, min)
#aggregate(Day.number ~ variable, all.data, max)
xout <- seq(342, 2651, by=3.35)
range(xout)

## [1]  342.00 2650.15

all.data1 <- na.omit(all.data)
mt <- plyr::dlply(all.data1,
                  "variable",
                  function(xx) pracma::interp1(x=xx$Day.number,
                                               y=xx$value,
                                               xi=xout,
                                               method="cubic"))
## Aside: the duplicated date that was previously fixed was ponly discovered by a warning message
## given by the pracma::interp1 function!!!

mt <- as.data.frame(mt)
mt <- cbind(Day.number=xout, mt)
mt <- gather(mt, variable, value, 2:13)
#ggplot(mt, aes(x=Day.number, y=value)) + facet_wrap(~variable, scales="free") + geom_line()

Fourth root transform

Next, because the original time series showed many sharp spikes, the time series were rescaled using a fourth-root power transformation (Fig. S1b). The sharp spikes bias “direct method” estimates of the Lyapunov exponent, because nearby pairs of reconstructed state vectors mostly occurred in the troughs between spikes. The average rate of subsequent trajectory divergence from these pairs is therefore an estimate of the local Lyapunov exponent in the troughs, which may be very different from the global Lyapunov exponent. By making spikes and troughs more nearly symmetric, the power transformation resulted in a much more even spread of nearby state vector pairs across the full range of the data for all functional groups in the food web. The transformation is also useful for fitting nonlinear models of the deterministic skeleton (used for nonlinear predictability and indirect method estimates of the Lyapunov exponent), which was done by least squares and therefore is most efficient when error variances are stabilized. Fourth-root transformation is intermediate between the square-root transformation that would approximately stabilize the measurement error variance in count data from random subsamples, and the log transformation that is usually recommended for stabilizing process noise variance due to stochastic variation in birth and death rates.

mt$fr.value <- mt$value^0.25

Detrend

The time series were then detrended using a Gaussian kernel with a bandwidth of 300 days (red line in Fig. S1b), to obtain stationary time series. Most species did not show long-term trends, except for the bacteria, detritivores (ostracods and harpacticoid copepods), dissolved inorganic nitrogen and soluble reactive phosphorus. One possible explanation for these trends in the microbial loop could be the slow accumulation of refractory organic material in the mesocosm, but we have not measured this component.

detrended <- group_by(mt, variable) %>%
  mutate(smoothed=ksmooth(Day.number, fr.value, kernel="normal", bandwidth=300)$y)
detrended$dt.value <- detrended$value - detrended$smoothed

Finally, the time series were linearly rescaled to have zero mean and a standard deviation of 1 (Fig. S1c).

final <- group_by(detrended, variable) %>%
  mutate(y=as.numeric(scale(dt.value)))
summarise(final, mean=mean(y), sd=sd(y))

## Source: local data frame [12 x 3]
## 
##                              variable          mean sd
## 1                          Cyclopoids -6.181296e-17  1
## 2                   Calanoid.copepods -3.437216e-18  1
## 3                            Rotifers -5.937068e-18  1
## 4                            Protozoa  2.108695e-18  1
## 5                   Nanophytoplankton  9.258537e-19  1
## 6                   Picophytoplankton -9.953683e-18  1
## 7                 Filamentous.diatoms  5.568714e-18  1
## 8                           Ostracods  1.266674e-17  1
## 9                       Harpacticoids -2.246616e-17  1
## 10                           Bacteria -2.384915e-17  1
## 11 Total.dissolved.inorganic.nitrogen -3.029117e-17  1
## 12        Soluble.reactive.phosphorus -3.780000e-17  1

We now have a y variable to work with in data.frame “final”, variable name “y”.

glimpse(final)

## Observations: 8280
## Variables:
## $ Day.number (dbl) 342.00, 345.35, 348.70, 352.05, 355.40, 358.75, 362...
## $ variable   (fctr) Cyclopoids, Cyclopoids, Cyclopoids, Cyclopoids, Cy...
## $ value      (dbl) 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e...
## $ fr.value   (dbl) 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.0...
## $ smoothed   (dbl) 0.3058059, 0.3092420, 0.3127101, 0.3161899, 0.31968...
## $ dt.value   (dbl) -0.3058059, -0.3092420, -0.3127101, -0.3161899, -0....
## $ y          (dbl) -0.4631805, -0.4842889, -0.5055944, -0.5269711, -0....

Zero removal

The time series of cyclopoid copepods, protozoa, filamentous diatoms, harpacticoid copepods and ostracods contained long sequences of zero values. This does not imply that these species were absent from the food web during these periods, but that their concentrations were below the detection limit. Time series dominated by many zeros can bias the statistical analysis. Therefore, these time series were shortened to remove long sequences of zero values, before the data transformation. The transformed data of all species in the food web are shown in Figure S2.

This is not done. May need to be done only for analyses for Table 1.

Figure S1 (visualising the transformation)

Choose a species to plot:

soi <- "Bacteria"

Raw and interpolated data:

g1 <- ggplot(filter(all.data, variable==soi), aes(x=Day.number, y=value)) +
  facet_wrap(~variable, ncol=2, scales="free_y") +
  geom_point(size=1, col="black") + geom_line(size=0.1) +
  scale_colour_manual(values = species.colour.mapping) + ggtitle("Raw and interpolated data")
g2 <- geom_line(data=filter(mt, variable==soi), aes(x=Day.number, y=value), size=0.25, col="blue")
g1 + g2

## Warning in loop_apply(n, do.ply): Removed 17 rows containing missing values
## (geom_point).

Fourth root transformed with trend:

g1 <- ggplot(filter(detrended, variable==soi), aes(x=Day.number, y=fr.value)) +
  facet_wrap(~variable, ncol=2, scales="free_y") +
  geom_point(size=0.5, col="black") + geom_line(size=0.1) +
  scale_colour_manual(values = species.colour.mapping) + ggtitle("Quarter power trans. and trend")
g2 <- geom_line(data=filter(detrended, variable==soi), aes(x=Day.number, y=smoothed), size=0.25, col="blue")
g1 + g2

Detrended and normalised:

g1 <- ggplot(filter(final, variable==soi), aes(x=Day.number, y=y)) +
  facet_wrap(~variable, ncol=2, scales="free_y") +
  geom_point(size=0.5, col="black") + geom_line(size=0.1) +
  scale_colour_manual(values = species.colour.mapping) + ggtitle("Detrended and normalised")
g1

Spectral analyses

Jason, Frank.

# Raw spectrum
spectra <- final %>% group_by(variable) %>% do(spectra = spectrum(ts(data=.$y, end=2650.15, deltat=3.35), log='no', method="pgram", detrend=F, plot=F))
spec <- spectra %>% do(data.frame(spec = .$spec[[2]], freq = .$spec[[1]], group = .[[1]]))

ggplot(spec, aes(y=spec, x=1/freq, group=group)) + geom_line() + facet_wrap(~group) +
coord_cartesian(ylim=c(0,40), xlim=c(0,240))

freq.est <- spec %>% group_by(group) %>% mutate(max_spec = max(spec), freq = freq)

## Warning: Grouping rowwise data frame strips rowwise nature

freq.est <- subset(freq.est, max_spec==spec, select=c(freq,group))
freq.est$freq <- 1/freq.est$freq
freq.est

## Source: local data frame [12 x 2]
## Groups: group
## 
##         freq                              group
## 1  2412.0000                         Cyclopoids
## 2   482.4000                  Calanoid.copepods
## 3   344.5714                           Rotifers
## 4    60.3000                           Protozoa
## 5   114.8571                  Nanophytoplankton
## 6  2412.0000                  Picophytoplankton
## 7  2412.0000                Filamentous.diatoms
## 8  1206.0000                          Ostracods
## 9  2412.0000                      Harpacticoids
## 10 2412.0000                           Bacteria
## 11 2412.0000 Total.dissolved.inorganic.nitrogen
## 12  301.5000        Soluble.reactive.phosphorus

# Welch's periodogram

wspectra <- final %>% group_by(variable) %>% do(spectra = pwelch(ts(data=.$y, end=2650.15, deltat=3.35), window=5, method="pgram", plot=F))
wspec <- wspectra %>% do(data.frame(spec = .$spec[[2]], freq = .$spec[[1]], group = .[[1]]))

ggplot(wspec, aes(y=spec, x=1/freq, group=group)) + geom_line() + facet_wrap(~group) +
coord_cartesian(ylim=c(0.1,100), xlim=c(0,240))+
scale_y_continuous(trans="log")

freq.est <- wspec %>% group_by(group) %>% mutate(max_spec = max(spec), freq = freq)

## Warning: Grouping rowwise data frame strips rowwise nature

freq.est <- subset(freq.est, max_spec==spec, select=c(freq,group))
freq.est$freq <- 1/freq.est$freq
frequency(final$y)

## [1] 1

ts <- as.ts(final$y, frequency = 0.3)
#time(ts)

Reproducing Table 1 using ELE supplement data.

Kevin, Marco. Use data.frame final, response variable y. Create dataset with zeros removed for this table:

final_nozeros <- final
final_nozeros$y <- ifelse(final_nozeros$y!=0, final_nozeros$y, NA)
# make it wide format
#library(tidyr)
final_long <- final_nozeros[final_nozeros$variable!="Cyclopoids" & final_nozeros$variable!="Filamentous.diatoms",c(1,2,7)]
final_wide <- spread(final_long, variable,y)
# we also removed data on Cyclopoids, not used in the table

Calculate correlation coefficients:

cor.coefs <- cor(final_wide[,c(-1)])

Only keep the upper triangle of the cor.pvals matrix:

for(i in 1:10){
  for(j in 1:10){
  cor.coefs[i,j] <- ifelse(i<j, cor.coefs[i,j], NA)
}}

Get p-vals too:

# https://stat.ethz.ch/pipermail/r-help/2005-July/076050.html
pn <- function(X){crossprod(!is.na(X))}
cor.prob <- function(X){
# Correlations Below Main Diagonal
# Significance Tests with Pairwise Deletion
# Above Main Diagonal
# Believe part of this came from Bill Venables
pair.SampSize <- pn(X)
above1 <- row(pair.SampSize) < col(pair.SampSize)
pair.df <- pair.SampSize[above1] - 2
R <- cor(X, use="pair")
above2 <- row(R) < col(R)
r2 <- R[above2]^2
Fstat <- (r2 * pair.df)/(1 - r2)
R[above2] <- 1 - pf(Fstat, 1, pair.df)
R
}
cor.pvals <- cor.prob(final_wide[,c(-1)])

Only keep the upper triangle of the cor.pvals matrix:

for(i in 1:10){
  for(j in 1:10){
  cor.pvals[i,j] <- ifelse(i<j, cor.pvals[i,j], NA)
}}

Add significance “stars” to cor.coefs from cor.pvals

cor.stars <- cor.pvals
cor.stars <- ifelse(cor.pvals<0.0001, "***",
                    ifelse(cor.pvals<0.001, "**",
                           ifelse(cor.pvals<0.05, "*", "")))

cor.cp <- cor.coefs
for(i in 1:10){
  for(j in 1:10){
  cor.cp[i,j] <- paste(round(cor.coefs[i,j],3), cor.stars[i,j])
}}

Remove NAs:

for(i in 1:10){
  for(j in 1:10){
  cor.cp[i,j] <- ifelse(cor.cp[i,j]=="NA NA", "", cor.cp[i,j])
}}
for(i in 1:10){
  for(j in 1:10){
  cor.cp[i,j] <- ifelse(i==j, "1", cor.cp[i,j])
}}

colnames(cor.cp) <- c("Calanoid.copepods  ","Rotifers  ","Protozoa  ","Nanophytoplankton  ","Picophytoplankton  ",
"Ostracods  ","Harpacticoid.copepods  ","Bacteria  ","Nitrogen  ","Phosphorus  ")
rownames(cor.cp)<-colnames(cor.cp)

Make it a table:

library(knitr)
table1b <- kable(cor.cp, format="html", col.names = colnames(cor.cp), align="c",
                caption="Table 1.'Correlations between the species in the food web. Table entries show the product–moment correlation coefficients, after transformation of the data to stationary time series (see Methods). Significance tests were corrected for multiple hypothesis testing by calculation of adjusted P values using the false discovery rate.' Significant correlations are indicated as follows: *: P<0.05; **: P<0.01; ***: P<0.001. 'The correlation between calanoid copepods and protozoa could not be calculated, because their time series did not overlap. Filamentous diatoms and cyclopoid copepods were not included in the correlation analysis, because their time series contained too many zeros.' (Beninca et al. 2008)")

table1b

Table 1.‘Correlations between the species in the food web. Table entries show the product–moment correlation coefficients, after transformation of the data to stationary time series (see Methods). Significance tests were corrected for multiple hypothesis testing by calculation of adjusted P values using the false discovery rate.’ Significant correlations are indicated as follows: : P<0.05; **: P<0.01;** : P<0.001. ‘The correlation between calanoid copepods and protozoa could not be calculated, because their time series did not overlap. Filamentous diatoms and cyclopoid copepods were not included in the correlation analysis, because their time series contained too many zeros.’ (Beninca et al. 2008)
	Calanoid.copepods	Rotifers	Protozoa	Nanophytoplankton	Picophytoplankton	Ostracods	Harpacticoid.copepods	Bacteria	Nitrogen	Phosphorus
Calanoid.copepods	1	-0.094 *	-0.032	-0.025	0.178 ***	0.004	-0.018	0.117 *	-0.041	-0.008
Rotifers		1	-0.006	-0.055	-0.037	0.163 ***	0.144 **	0.282 ***	-0.046	0.147 **
Protozoa			1	0.044	-0.018	-0.023	-0.025	-0.025	-0.005	-0.005
Nanophytoplankton				1	-0.033	-0.024	-0.036	-0.01	-0.045	-0.047
Picophytoplankton					1	-0.018	-0.054	0.033	-0.034	-0.029
Ostracods						1	0.289 ***	0.211 ***	-0.12 *	0.2 ***
Harpacticoid.copepods							1	0.293 ***	-0.035	0.266 ***
Bacteria								1	-0.129 **	0.343 ***
Nitrogen									1	-0.071
Phosphorus										1

# differs from the one published by Beninca et al.!

Predictability (Figure 2)

To be done.

Use data.frame final, response variable y.

Lyapunov exponents by direct method (Figure 3)

Gian Marco, Mikael.

Use data.frame final, response variable y.

Lyapunov exponents by indirect method

Dennis, Vanessa.

Use data.frame final, response variable y.

Reproduce anaylses of Beninca et al (2008)

Owen Petchey

29 Jan 2015

Introduction

First get the raw data into R and tidy it.

Reproducing figure 1b through 1g

Data transformation

Interpolation

Fourth root transform

Detrend

Zero removal

Figure S1 (visualising the transformation)

Spectral analyses

Reproducing Table 1 using ELE supplement data.

Predictability (Figure 2)

Lyapunov exponents by direct method (Figure 3)

Lyapunov exponents by indirect method